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Abstract. Recent studies of the fragmentation functions using the Nambu-Jona-Lasinio (NJL) - Jet 
model have been successful in describing the quark fragmentation functions to pions and kaons. The 
NJL-Jet model employs the integral equation approach to solve for the fragmentation functions in 
quark-cascade description of the hadron emission process, where one assumes that the initial quark 
has infinite momentum and emits an infinite number of hadrons. Here we introduce a Monte Carlo 
(MC) simulation method to solve for the fragmentation functions,, that allows us to relax the above 
mentioned approximations. We demonstrate that the results of MC simulations closely reproduce 
the solutions of the integral equations in the limit where a large number of hadrons are emitted 
in the quark cascade. The MC approach provides a strong foundation for the further development 
of the NJL-Jet model that might include many more hadronic emission channels with decays of 
the possible produced resonances, as well as inclusion of the transverse momentum dependence 
(TMD), all of which are of considerable importance to the experimental studies of the transverse 
structure of hadrons. 
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INTRODUCTION 

The novel efforts to extract the quark fragmentation functions from various experimental 
data [1,2] have generated a renewed interest in the long studied subject of analyzing hard 
scattering reactions [3, 4, 5, 6, 7, 8, 9]. The analysis of the transversity quark distribution 
functions [8, 10] and a variety of other semi-inclusive processes [11, 12] also critically 
depend on the knowledge of the fragmentation functions. 

The recent development of the NJL-Jet model [13, 14] has led to an ever more 
sophisticated model for calculating quark fragmentation functions in an effective chiral 
quark theory. Here the quark fragmentation is modeled as a quark-cascade process, 
depicted in Fig. 1, and the multiplicative ansatz of Ref. [3] was used to derive a set of 
coupled integral equations for the fragmentation functions. The advantages of the NJL- 
Jet model are the absence of ad-hoc parameters included, i.e. all the NJL parameters 
are fixed independent of any experimental data on fragmentation functions. Moreover, 
the momentum and isospin sum rules for the solutions of fragmentation functions are 
naturally satisfied. On the other hand, the model assumes that the initial quark has infinite 
momentum (Bjorken limit) and emits an infinite number of hadrons, so the probabilities 
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FIGURE 1. Quark cascade. 



scale with the fraction of the light-cone momentum left in the quark cascade, while in 
medium-energy experiments only a few hadrons are emitted per struck quark. Another 
problem appears as more and more possible emission channels are included in model, as 
the numerical task of solving the corresponding integral equations becomes challenging. 
Lastly, the inclusion of the transverse momentum dependence of the fragmentation 
functions seems unachievable in the integral equation approach. 

Here we propose to use Monte Carlo method to solve for the fragmentation functions 
within the quark-cascade model. Our goal is to demonstrate the viability of the approach 
of replacing the integral equations by providing very similar order of precision in deter- 
mining the fragmentation functions over the entire range of the hght-cone momentum 
fraction z of the initial quark carried by the emitted hadrons. 

MONTE-CARLO APPROACH AND THE ASYMPTOTIC 
EQUIVALENCE TO THE INTEGRAL EQUATION METHOD 

The NJL-Jet model ([13, 14]) relied on a set of coupled integral equations to solve for 
the fragmentation functions: 

J^\{z)dz = d\{z)dz + Y. CdQ{y)dyD\ f^) ^, (1) 
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where we assumed that the quark has infinite momentum and produces an infinite 
number of hadrons. Here J^\{z) denotes the fragmentation function of quark q to hadron 

h carrying light-cone momentum fraction z, dqiz) is the "elementary" fragmentation 
function of quark q emitting hadron h at each vertex in the quark-cascade and the sum 
on the right hand side is over all active flavors of quarks included in the model. 

Here we propose to calculate the fragmentation functions using Monte-Carlo (MC) 
simulations akin to the method described in the Ref. [15] using the probabilistic inter- 
pretation of the fragmentation functions: Dq{z) is the probability to emit a hadron h 
carrying the light-cone momentum fraction z to z + dz of initial quark ^ in a quark-jet 
picture. The quark goes through a cascade of hadron emissions, where at every emission 
vertex we choose the type of emitted hadron h and its light-cone momentum fraction z 
(of the fragmenting quark) by randomly sampling the corresponding probabilities of the 
elementary fragmentations, dq{z) that are calculated within the NJL model (in general 
these can be calculated in any effective quark model). We keep track of the flavor and 
the light-cone momentum fraction of the initial quark left to the remaining quark, also 
recording the type and light-cone momentum fraction of the initial quark transferred to 



the emitted hadron. We stop the fragmentation chain after the quark has emitted a prede- 
fined number of hadrons, Nunks (Alternatively, we can stop the chain after the remnant 
quark in the cascade has less than a given fraction of the initial quark's light-cone mo- 
mentum, ZMin)- We repeat the calculation Ngims times with the same initial quark flavor, 
q until we have sufficient statistics for the emitted hadrons. We extract the fragmenta- 
tion functions by calculating the average number of hadrons of type h with light-cone 
momentum fraction z to z-|-Az, (A/^'(z,z-f Az)) and expressing them in terms of frag- 
mentation functions: 
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From the construction it is obvious that the fragmentation functions calculated using 
the integral Eq. (1) should be equivalent to those calculated using the MC method in 
the limit Nunks °° and Nsims °°- The plots in Fig. 2 show that the solutions for 
fragmentation functions from both methods are indeed equivalent with a large enough 
number of emitted hadrons within statistical errors. 
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FIGURE 2. Comparison of the solutions for quark fragmentation function Df^ (z) in NJL-Jet model 
with only non-strange pseudoscalar mesons calculated from integral equations Eq. (1) and MC simulation. 

MC also allows us to study the dependence of the resulting fragmentation functions on 
the number of hadrons emitted by the quark in the cascade, which could well be relevant 
to many medium-energy experiments. The plot in Fig. 3 show that the solution for 
zD^^ (z) with N links = 1 (equivalent to the elementary fragmentation function zd^^ (z)) 
is peaked at z ~ 0.8. As we increase the number of emitted hadrons, the solution 
increases in the low z region due to the hadrons emitted further in the quark jet, where 



the fragmenting quark typically has a small fraction of the initial quark's light-cone 
momentum. We can readily see that the solutions saturate after including only a few 
emitted hadrons, where there is virtually no difference between solutions with Nunks — 8 
and Nunks = 20, and the discrepancy with the solution of the integral equations only 
occurs at vanishingly small values of z- Thus we can reliably use the solutions of MC 
simulations with Nunks > 8. 
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FIGURE 3. The dependence of the solutions for zD^^ on Nunks- 

FINAL REMARKS 

In the current article we demonstrated that the Monte Carlo approach to calculating 
the fragmentation functions in NIL- Jet framework is a powerful and reliable method. 
We reproduced the fragmentation functions calculated as solutions of the previously 
employed integral equations, where only the light quarks and pions were included. 
Moreover, we showed that the MC approach allows for the flexibility to surpass the 
model limitations necessary in formulating the integral equations, i.e. in the future MC 
studies we can assume an initial quark carrying only a finite momentum, thus emitting a 
finite number of hadrons. 

A further advantage of the MC approach is in reducing the numerical task in solv- 
ing for the fragmentation functions when including many more channels for emitted 
hadrons. Here solving the integral equations requires inverting larger and larger matrix, 
while the MC procedure can be drastically sped up by trivially parallelizing the task and 
solving simultaneously on computer clusters. 

It is clear that for a more complete description of the quark fragmentation both 
vector meson and nucleon anti-nucleon channels, as well as the strong decays of the 




produced resonances need to be included in the calculations. This is because as the high 
z region of the fragmentation functions are dominated by "few-step" transitions where 
the availability of the additional fragmentation channels might have a noticeable effect. 
The future development of the NJL-Jet model would also allow access to the transverse 
momentum distribution of the produced hadrons, thus becoming relevant for the analysis 
of a large variety of semi-inclusive data. The MC approach provides a strong platform 
for implementing these and other possible extensions of the NJL-Jet model that would 
allow for a much more detailed description of the physical picture. 
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